function [nFs, mXY, vShut, bReload] = Frequency_Adaptation(bAutoLoad)
% GalvoScanControl
% Scan script template: Frequency_Adaptation
%
% Description:
% Run many pulse trains, at different frequencies to test frequency and
% adaptation properties. The order of presented frequencies is (by default)
% randomized.
%
% Parameters:
%   train duration (pulses)
%   train repetitions (#)
%   pulse frequencies (Hz)
%   pulse duration (ms)
%   voltage (V)
%

% Scan scripts should return four variables (in this order):
%  nFs      The sample rate (1/sec between scan points)
%  mXY      A matrix with two columns in the format [X Y], where
%               X  are x-positions
%               Y  are y-positions
%           The unit of this matrix is millimeters, relative to the origin
%  vShut    Vector with same length as mXY where values indicate shutter
%           open/close state (when a shutter is controlled) or the beam
%           intensity (when the output is directed to an intensity
%           modulated laser module, e.g. Coherent CUBE).
%           Note: Values of vShut should not exceed +/- 5V. Values beyond
%           these limits will be fixed at +/- 5V by the software.
%  bReload  Boolean indicating whether the script should be re-loaded each
%           time before executing. Setting this to true (1) will decrease
%           the repetition rate and should only be used for long programs
%           that are not executed at high frequency.

persistent p_sCycles p_sTrainReps p_sFreqs p_sPulseDur p_sVolts
bReload = 1;

% Default values
if isempty(p_sCycles), p_sCycles = '10'; end
if isempty(p_sTrainReps), p_sTrainReps = '6'; end
if isempty(p_sFreqs), p_sFreqs = '[2 5 10 20 30 40 50]'; end
if isempty(p_sPulseDur), p_sPulseDur = '2'; end
if isempty(p_sVolts), p_sVolts = '5'; end

if ~bAutoLoad
    cAnswers = inputdlg({'Train duration (pulses)', ...
        'Train repetitions (#)', ...
        'Train frequencies (Hz)', ...
        'Pulse duration (ms)', ...
        'Intensity (0-5 V)'}, ...
        'Frequency Adaptation', ones(1, 5), ...
        {p_sCycles, p_sTrainReps, p_sFreqs, p_sPulseDur, p_sVolts});
    if isempty(cAnswers), return, end
    
    % For some reason, closing the input dialog takes forever.
    % Placing this pause here prevents this.
    pause(0.1)
    
    p_sCycles = cAnswers{1};
    p_sTrainReps = cAnswers{2};
    p_sFreqs = cAnswers{3};
    p_sPulseDur = cAnswers{4};
    p_sVolts = cAnswers{5};
end

nCycles = str2double(p_sCycles);
nTrainReps = str2double(p_sTrainReps);
vFreqs = eval(p_sFreqs);
nPulseDur = str2double(p_sPulseDur);
nVolts = str2double(p_sVolts);

% Truncate intensity to 0 - 5 V
if nVolts < 0, nVolts = 0;
elseif nVolts > 5, nVolts = 5; end

% Duration between trains is 1 / lowest frequency
nTrainSep = 1/min(vFreqs);

% Increase frequency vector by the number of train repetitions
vFreqs = repmat(vFreqs, 1, nTrainReps);

% Randomize order of frequency presentation
vFreqs = vFreqs(randperm(length(vFreqs)));

% Create a pulse train for each frequency
vShut = [];
for f = 1:length(vFreqs)
    % Generate analog trace from pulse train parameters
    nFreq = vFreqs(f);
    nPeriod = (1 / nFreq) * 1000; % ms
    nPulseDown = nPeriod - nPulseDur;
    
    % Single pulse
    vOpen = [ones(nPulseDur, 1)];       % Open period
    vClosed = [zeros(nPulseDown, 1)];   % Close period (0V)
    
    % Compose pulse train
    vShut_this = [];
    for i = 1:nCycles
        vShut_this = [vShut_this; vOpen .* nVolts; vClosed];
    end
    
    % Insert into global analog vector
    vShut = [vShut; vShut_this];
    
    % Insert spacer between trains (1 / lowest frequency)
    if f < length(vFreqs) % don't insert after last trial
        vShut = [vShut; zeros(nTrainSep * 1000, 1)];
    else
        vShut = [vShut; 0];
    end
    
end

% Position of beam is constant throughout this protocol
mXY = repmat(0, length(vShut), 2);
nFs = 1000;

return

